Symbolic dynamics of biological feedback networks 
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We formulate general rules for a coarse-graining of the dynamics, which we term 'symbolic dynam- 
ics', of feedback networks with monotone interactions, such as most biological modules. Networks 
which are more complex than simple cyclic structures can exhibit multiple different symbolic dy- 
namics. Nevertheless, we show several examples where the symbolic dynamics is dominated by a 
single pattern that is very robust to changes in parameters and is consistent with the dynamics being 
dictated by a single feedback loop. Our analysis provides a method for extracting these dominant 
loops from short time series, even if they only show transient trajectories. 
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Many biological systems can be described by directed 
networks, where nodes represent different components 
and arrows represent interactions. In cell biology, nodes 
are molecules, while arrows stand for complex forma- 
tion, protein modification, transcription regulation, etc. 
Ecosystems constitute another example, where nodes are 
species and arrows represent predation, competition and 
symbiosis. Biological functions are often performed by 
specific small subnetworks, or modules [1]. A dynamical 
model of a module requires, beyond the knowledge of the 
network structure, some hypothesis on the form of the 
interactions, which are often poorly characterized. It is 
then crucial to develop techniques to study the qualita- 
tive dynamics of modules assuming limited information. 

We present a method to obtain information about pat- 
terns in the dynamics of biological feedback networks. 
Given the network structure, i.e., which nodes acti- 
vate/repress which other nodes, it is possible to pre- 
dict the ordering of maxima and minima of the dynam- 
ical variables. Vice versa, from an experimentally ob- 
served ordering one can obtain some information about 
the structure of the network. The method is a general- 
ization of the one introduced in [2| for the dynamics of 
a single negative feedback loop, without any cross links, 
where a unique pattern is allowed. A More complex net- 
work structure [H, HI, H, S 0, S] allows multiple dynamical 
patterns. Moreover, a particular observed pattern could 
originate from different network structures. Neverthe- 
less, our method can reduce the possibilites and provides 
non-trivial information that can guide experiments. Our 
technique also applies when the dynamics is transient, so 
that information can be obtained, for example, by watch- 
ing how the concentrations of proteins and genes belong- 
ing to a given module relax to a stationary state after a 
perturbation. This extends the use of our formalism to 
cases in which oscillations are damped 0- 

We consider a system described by A dynamical vari- 
ables, Xi, i — 1 ... A, which we call "densities" and could 
represent, for example, the concentrations of the chemical 
species composing the network/module. We assume that 



they evolve with time in a deterministic way, according 
to a system of differential equations: 
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Many possible dynamical systems may correspond to 
a given network. The only constraints we impose are 
that the interactions be monotonic, i.e., each off-diagonal 
element of the Jacobian matrix, d x . gi , is either positive 
everywhere in phase space (when node j activates i), or 
negative everywhere (J represses i), or zero everywhere 
(no arrow from j to i). In words, activators are always 
activators and repressors are always repressors. Indeed, 
transcription factors rarely switch from being activators 
of a particular gene to repressors at different densities; 
a predator of a particular species does not become its 
prey when their abundances change. We do not require 
monotone self-interactions: a variable may activate or 
repress itself depending on the densities. 

We associate to each state (x±,X2 ■ ■ - Xn) a symbol 
such as (+,—,—,...+). This A-component sign vec- 
tor describes which densities are increasing and which 
are decreasing at a given time: the i-th component is 
just the sign of g%(x\, x% . . . Xn)- Such a representa- 
tion divides the phase space into sectors, each associ- 
ated with a symbol, in which each density has a defi- 
nite increasing/decreasing behavior. The sectors' bound- 
aries are the nullclines, i.e. the manifolds satisfying 
gi(%i, x% . . . xn) = 0. Our goal is to determine the con- 
ditions under which the trajectory can cross a nullclinc. 
This requires a density to change from increasing to de- 
creasing (or vice versa) and is equivalent to determining 
when the density can have a maximum or minimum. 

For example, a minimum for the variable Xi corre- 
sponds to a crossing of the nullcline gi = from the 
region gi < to the region gi > 0. This is possible 
only if, somewhere on the nullclinc, the scalar product 
between the vector field g and the vector V<?i (which is 
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FIG. 1: Simple example of the use of symbolic dynamics, 
(a) (left) Scheme of the network. We represent activation by 
a normal arrow, and repression by a barred arrow, (right) 
Corresponding transition network; we removed for clarity the 
symbols (+-+) and (-+-) which have no incoming links, (b) 
Dynamics of the 3 variables as a function of time with a weak 
cross-link (a = 10, see text), showing the transition cycle in 
solid arrows in (a). Inset shows the same on a log-scale, (c) 
Dynamics for a stronger cross-link (a — 50, see text) which 
includes the transitions shown by dashed arrows, zoomed in 
the inset. In all plots x\ is blue, X2 is green and X3 is red. 



normal to the nullclinc gi = 0) is positive: 



■%n) d Xj gi(xi,X2,...x N ) > 0. (2) 



The i = j term is excluded since it is zero on the 
nullcline. By assumption, all the derivatives have fixed 
signs, and in any given sector the gj's also have fixed 
signs (encoded in the associated symbol) . If the symbol 
and derivative signs are such that each term is negative, 
then the sum cannot be positive. This implies the rule: 
A variable cannot have a minimum if all its repressors 
are increasing and all its activators are decreasing. 

Similarly, for maxima: 
A variable cannot have a maximum if all its repressors 
are decreasing and all its activators are increasing. 

Using the above two rules we can construct a network 
of allowed transitions for a given biological module, with 
one node for each symbol and an arrow for each transi- 
tion that does not violate the above rules. Note that the 
transition network will only have arrows connecting sym- 
bols differing by a single sign, because each maximum or 
minimum corresponds to a single sign flip. 

We first consider the example network of Fig. la(left): 
a three-species negative feedback loop with a cross-link 
from node 3 to 2 that introduces a positive feedback. By 
checking all the allowed transitions we construct the cor- 



responding transition network, shown in Fig. la(right). 

For example, from the symbol ( h +) the transition 

( — h+) — > ( 1 ) is ruled out because all the acti- 
vators of node 3 are increasing, therefore it cannot have 
a maximum. Similarly we rule out all transitions from 

it except ( — h +) — > ( h). The result, in this case, 

is a simple modification of the transition network for a 
single negative feedback loop shown by the solid arrows 
111 Fig. la(right)@. With the cross link present, the 
additionally allowed transitions are the ones shown with 
dashed arrows. The following dynamical system illus- 
trates these possibilities: x\ = c — x\X^/(ki + xi); Xi = 
x\ + a [6(x3 — k 2 ) — 1] — x%; X3 = x 2 — x%. The major 
nonlinearity is the Heaviside s tep function: 6{x) = for 
x < 0, and 6(x) = 1 for x > 0[2fJ. By choosing parame- 
ters such that the cross link is weak (c = 30, a — 10, k\ — 
0.1, k 2 = 20) one obtains dynamics of Fig. lb, which is 
identical to the simple 3-species loop. As the strength of 
the cross link is increased (a = 50), the symbolic dynam- 
ics changes to also exhibit the dashed transitions. This is 
shown in Fig. lc where variable x 2 develops a new small 
maximum, thus changing the symbolic dynamics. 

We now move on to a system that exhibits a richer 
range of dynamical behaviours (see Fig. [2^). It con- 
sists of two negative feedback loops, coupled via a shared 
species. This network has been widely studied in the eco- 
logical literature 0, [13, II | as a model for three trophic 
level ecosystems: species X3 feeds on x 2 , and x 2 feeds 
on x\. The chaotic properties of this motif have been 
used to interpret data from the Canadian lynx-hare cy- 
cle, showing irregular oscillations 

We consider first the Hastings-Powell (HP) model [l(| 
as a dynamical system corresponding to this network: 

x\ = rx\ (1 — kx\) — o.\X\X 2 j{\ + b\Xi) 

x 2 = -d x x 2 + a-iX-iXz/il + bi^i) - a 2 x 2 x 3 /(l + b 2 x 2 ) 

x 3 = -^22:3 + a 2 x 2 x 3 /(l + b 2 x 2 ) (3) 

with the following parameter choices: ot\ — a 2 = 4, 
b\ = b 2 = 3, d\ = .4, d 2 = .6, k = 1.5. By increasing 
the parameter r, a stable limit cycle undergoes a series 
of period doubling bifurcations, followed by the onset of 
chaos. A projection of the attractor on the x 2 — x 3 plane 
is shown in Fig. ([3]). The chaotic trajectory looks similar 
to the periodic one, except for the irregular behavior of 
the amplitude llj. This means that the same sequence 



corresponding to the periodic orbit is observed after the 
onset of chaos. By increasing r even more, we found a 
regular window with a change in the symbolic dynamics 
(the "kick" , shown in red in the attractor in Fig. [3] and 
in the bifurcation diagram, Fig. |4^, and corresponding 
to the blue dashed transition in Fig. 2b). The kick is 
still present when, by further increasing r, the dynamics 
becomes chaotic again. 

The conclusion is that the same symbolic dynamics 
observed close to the Hopf bifurcation is found in a large 
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FIG. 2: Network of two coupled two-species oscillators, (a) 
Structure of the network, (b) The transition network for this 
3-node system. Black arrows indicate all the allowed transi- 
tions. Blue arrows are the transitions actually observed in the 
HP system and red arrows are the transitions observed in the 
BHS model (see text). Dashed arrows indicate "kicks", i.e., 
transitions which are not observed close to the Hopf bifurca- 
tion. 
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FIG. 3: Two dimensional projection of the attractor of the 
system of equation ((3} for different values of the control pa- 
rameter r = 2.0 (top-left), r — 2.6 (top-right), r — 3.0 
(bottom- left) r = 3.3 (bottom-right). 

region of parameter space. We compare the results with 
a different system corresponding to the same network, 
the model by Blausius, Huppert and Stone (BHS)[9(: 

x\ = x\ — a.\X\x%J{\. + kx%) 

x 2 = -dx 2 + atixix 2 /(l + kxi) - a 2 x 2 x 3 

x 3 = c(xl - x 3 ) + a 2 x 2 x 3 (4) 



with parameters ax = 2, a.^ = d = 1, k = 0.12, 
x 3 = 0.006. Here, a convenient control parameter is c. 
We observe the same scenario in the bifurcation diagram 
(see Fig. (0b)): periodic orbit, then chaotic but same 
periodic symbolic dynamics, then different symbolic dy- 
namics in a regular window and, finally, chaotic symbolic 
dynamics. Note, however, that the periodic symbolic dy- 
namics observed close to the Hopf bifurcation is different 
from that observed in the HP model. 

To test the robustness of the two sequences, we tried to 
change the functional form of the interaction between x 2 
and £3 by setting b 2 — in the HP model or, conversely, 
introducing saturated response in the BHS model. We 
also tried to vary the parameters of both systems, by 
up to 50% from their default values. The two symbolic 
sequences were not affected by any of these changes. A 
possible cause for this robust difference could be the logis- 
tic term in the first equation of ([3]) , acting as a regulator 
so that the full dynamics is bottom-up controlled. 

The difference between the symbolic dynamics of the 
HP and BHS systems can be used for model selection: in 
the example of the Canadian lynx system, one has access 
only to the lynx population time series, but temporal 
measurements of the hare and grass abundances could 
be used to understand which model is more appropriate. 
Interestingly, from the point of view of maxima/minima 
order, these two systems behave like two different, single 
negative feedback loops 0: 3 H 2 H H 3 (HP system) 
and 1 -> 2 -> 3 H 1 (BHS system). Both these "effective" 
loops would include an "effective" interaction between 
variables x\ and x%. 




FIG. 4: Bifurcation diagrams, (top) HP model (|3}, minima 
of X2 plotted versus r. (bottom) BHS model Q, maxima of 
x\ versus c. In both plots, red dots indicate the appearance 
of "kicks" in the trajectory and symbolic dynamics (see text). 

So far we have gone from a known network to the 
transition network to the time series. The reverse pro- 
cess uses the transitions observed in an experimental 
time series to infer information about the underlying net- 
work. For example, the circadian oscillations of the three 
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genes kaiA, B, C in cyanobacteria 13[ were shown in ref. 
[2| to have the following symbolic dynamics (B,A, C): 

(+ + +) -(- + +)-(--+)-( )_>(+__)_ 

(+H — ) — » (+ + +). Several networks are consistent with 
this pattern - the simplest is the loop B — > A — > C H £?, 
as suggested in ref. [2J. Experiments have shown that 
A — » C and CHS, and that all three genes are es- 
sential for oscillations 14} . With this information we can 
get non-trivial guidelines for which further interactions to 
look for experimentally: (i) kaiA must be either activated 
by kaiB, or repressed by kaiC (or both), (ii) if kaiA is 
not activated by kaiB then, in addition to kaiC H kaiA, 
kaiB must activate kaiC, so that the underlying network 
looks similar to Fig. 2a. Of course, these predictions are 
for "effective" interactions, which at the molecular level 
could involve multiple intermediates, such as chemical 
complexes and various protein activity states. Ref. Q 
shows how the method can reconstruct effective interac- 
tions even in the presence of intermediate species. 

This circadian example also points out how much infor- 
mation our method provides. The transition (+ + +)—> 
( — h +) means that either B H A or C H A. Later tran- 
sitions show that either A — > B or C H B, and A — ► C 
or B — > C . Even without the extra experimental in- 
formation, our method reduces the possibilities for the 
adjacency matrix of the underlying network from 3 6 to 
5 3 , a factor of ks 6. In a general N node system, with 
M independent observed transitions, the fraction of al- 
lowed adjacency matrices is [1 — (2/3) ^J ; the smaller 
the network and the more the transitions seen, the more 
useful the method. A full oscillation cycle would show at 
least N independent transitions. If the system instead 
reaches a fixed point, the transient can still be used. 

Our method can be considered as complementary to 
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the "threshold method", described in Refs 
[Tij l , which provides a different way of dividing the phase 
space into sectors, based on a choice of thresholds for 
each variable. The "threshold" method generates a tran- 
sition diagram, which depends on parameter values. It is 
particularly suited to cases where the input functions are 
Boolean or step-like, so thresholds can be easily identi- 
fied, and self-interactions are piecewise linear [19']. Our 
"derivative" method, in contrast, requires no choice of 
thresholds, generates a diagram independent of parame- 
ter values, and works for arbitrary self-interactions, but 
requires monotonicity in the other interactions. 

In summary, we showed a method to construct a sym- 
bolic transition network that imposes a strong constraint 
on the dynamics monotone systems, like many biological 
modules. In all the cases we studied, the periodic sym- 



bolic dynamics observed close to the Hopf bifurcation 
is found in a large region of parameter space, even when 
the system becomes chaotic. This explains the commonly 
observed phenomenology of a chaotic attractor consist- 
ing of oscillations with randomly varying amplitude [H| . 
The oscillatory systems we looked at produce a symbolic 
sequence identical to that of a single negative feedback 
loop for most studied parameter values. By identifying 
these loops, our method can be used to derive minimal 
models of complex, oscillatory biological systems. 

We acknowledge Leon Glass for stimulating discussions 
and the Danish National Research Foundation and VIL- 
LUM KANN RASMUSSEN Foundation for funding. 
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We can safely introduce a discontinuous field because our 
argument works as long as trajectories are continuous. 



